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Abstract. We propose and investigate a numerical shooting method for computing geodesies in 
the Weil-Petersson (WP) metric on the universal Teichmiiller space T(l). This space, or rather 
the coset subspace PSL2(lR)\Diff(5 1 ), has another realization as the space of smooth, simple 
closed planar curves modulo translations and scalings. This alternate identification of T(l) is 
a convenient metrization of the space of shapes and provides an immediate application for our 
algorithm in computer vision. The geodesic equation on T(l) with the WP metric is EPDiff(5' 1 ), 
the Euler-Poincare equation on the group of diffeomorphisms of the circle S 1 , and admits a class of 
soliton-likc solutions named Teichons |20| . Our method relies on approximating the geodesic with 
these teichon solutions, which have momenta given by a finite linear combination of delta functions. 
The geodesic equation for this simpler set of solutions is more tractable from the numerical point 
of view. With a robust numerical integration of this equation, we formulate a shooting method 
utilizing a cross-ratio matching term. Several examples of geodesies in the space of shapes are 
demonstrated. 
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1. Introduction 

We consider the Weil-Petersson (WP) metric on the coset space PSL2(M)\DifT(S' 1 ). This coset 
space (or its completion in the WP metric or in the Teichmiiller topology) is known as the universal 
Teichmiiller space and is well-known in many contexts: in the classification of Riemann surfaces [17] . 
conformal and quasi-conformal maps [23) . string theory [5] and most recently computer vision [30] . 
Its completion in the WP metric is an infinite dimensional homogeneous complex Kahler-Hilbert 
manifold [31] ■ 

As we will explain in Section [2] below, a particular dense subset of the universal Teichmiiller 
space T(l) is given by PSL 2 (IR)\Diff (S 1 ), where Diff (S* 1 ) is the group of C°° diffeomorphisms of 
S 1 , and PSL^M) is a subgroup of the Mobius selfmaps of the unit disk, see Q and the surrounding 
discussion. This coset space is a Riemannian manifold for the WP metric and has another realization 
as the space of smooth simple closed curves modulo translations and scalings. (Therefore we will 
hereafter use the terms 'shape', 'diffeomorphism', 'fingerprint', or 'welding map' to refer to members 
in this dense subset of T(l).) Endowing a shape space with a metric and computing geodesies 
between shapes is a central problem in computer vision. It aids in recognition and classification, 
enables computational strategies to address the clique problem, and allows us to perform statistics on 
shapes. Another application is in the emerging field of computational anatomy, where quantitative 
analysis of anatomical variability is important |14j . It is this major application that provides the 
motivation for this work. 

There are several advantages to using the Weil-Petersson metric to compare two-dimensional 
shapes. First, any two smooth shapes can be connected with a Weil-Petersson geodesic [13] ■ Second, 
all sectional curvatures of the metric are negative [3T]. Thus geodesies connecting two shapes are 
unique (18) . We are not aware of any other metric currently used in the pattern theory literature 
that has these properties. Ultimately, the uniqueness of geodesies allows one to perform consistent 
statistical analysis on shape databases via the initial momentum representation of the shape, but 
this application is beyond the scope of this paper. 

Geodesic equations of groups of diffeomorphisms on a manifold M were first studied in Arnold's 
ground-breaking paper |3] . Arnold considered in particular the group of volume preserving diffeo- 
morphisms of Euclidean space in its L 2 metric and found the geodesic equation for the vector field 
v(x,t) to be Euler's fluid flow equation (see [4] for a full exposition). Other examples include 
the periodic Korteweg-deVries (KdV) equation and the periodic Camassa-Holm (C-H) equation [7J. 
These equations are geodesic equations on the Virasoro group, a central extension by S 1 of the 
group Diff(S' 1 ) of the diffeomorphisms of S 1 , for the L 2 and H 1 metric respectively. KdV and 
C-H are two completely integrable partial differential equations and have soliton solutions. Holm 
and collaborators have found that the geodesic equation on Diff(R") admits special solutions with 
many of the properties of solitons: for each fixed time, they are diffeomorphisms which are largely 
localized in space and retain their general shape as they evolve; furthermore they interact somewhat 
like KdV solitons )12) . There are not, however, infinitely many conserved quantities so they are not 
true solitons. 

Singular solutions first arose as peakons (from 'peaked solitons') for a completely integrable 
Hamiltonian water wave equation, C-H in jJJ . The peaks occurred where the velocity profiles of the C- 
H equation had discontinuity in its slope. These peaks correspond to Dirac delta distributions of the 
associated momentum. The EPDiff equation for other metrics was later found independently in [35] , 
and its singular solutions were shown to be important as landmarks in shape analysis |141 127] . Later 
they were shown to comprise a singular momentum map for the right action of the diffeomorphisms 
on embeddings in any dimension [15]. Currently, the use of EPDiff and its landmark solutions is 
standard in shape analysis [T^l [H] [55] . 
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It turns out that considering the Weil-Petersson metric on the coset space PSL2(R)\Diff (S 1 ) 
yields another example of a geodesic equation that is similar to KdV and C-H. This equation 
describing evolution of the velocity field v(t, 8) is 

(1) m t + 2mve + vmg — 0, where m = —H(vg + vgee), 

and T-L is the periodic Hilbert transform defined by convolution with 7^ctn(#/2). 

It is not known if (JlJ is completely integrable but it admits a class of soliton-like solutions which 
we consider in this paper: solutions in which m can be represented as a finite sum of weighted 
Dirac delta functions. Darryl Holm suggested the portmanteau teichons to describe these soliton- 
like solutions on Teichmuller space and their corresponding geodesies. We adopt this terminology 
in this paper. 

We use an iV-teichon ansatz (a sum of N teichons) to reduce the integro-diffcrcntial equation 
([I]) to a finite-dimensional system of ordinary differential equations. In this way we approximate 
geodesies between any two points in PSL 2 (IR)\Diff (S* 1 ) with 7V-teichon geodesic evolutions. We 
use this teichon formulation to shoot from an initial shape to a terminal shape that must then be 
compared with the target shape. Because we are considering the coset space PSL2(R)\Diff (S 1 ), 
using a standard matching term on diffeomorphisms is not possible: the quotient space ambiguity 
prevents straightforward comparisons of diffeomorphisms (e.g. pointwise matching). We address 
this difficulty by using cross-ratios in the matching term; their invariance within an equivalence 
class on PSL2(M)\Diff (S 1 ) allows for accurate matching on the coset space. 

One existing geodesic computing algorithm, described in [30] and investigated in [21], suffers 
from several limitations. In particular there are numerical difficulties in matching shapes that are 
not close to a circular shape. The current approach aims to mitigate this shortcoming. A recent 
approach given in [11] is more competitive with our approach. 

This paper is organized as follows. Section [2] introduces the background on universal Teichmuller 
space T(l), fingerprints (also called welding maps), and the Weil-Petersson metric. With the WP 
metric, Section [3] discusses the geodesic equation on the Teichmuller space (also known as EPDiff), 
and Section [4] discusses teichon solutions of EPDiff. In Section [5] we describe the shooting method 
and the matching functional used for shooting, along with details of the gradient computation. In 
Section [6] we demonstrate the utility of our method with several examples. 

2. Shapes as diffeomorphisms of the circle S 1 

2.1. Fingerprints. Let Di nt be the open unit disk in the complex plane C, i.e. Dj nt = {z £ C | \z\ < 
1}, and let D ext = {z e C | \z\ > 1} be its exterior. For every simple closed curve T in C denote by 
Tint its union with the region enclosed by it, and denote by r ext its union with the infinite region 
outside of T (including oo). 

Then by the Riemann mapping theorem, for all T there exist two conformal maps 

/int : D int — > r int , 
/ext '■ D oxt — > Text- 

The interior map /; n t is unique up to replacing f- mt by /i n t o A for any Mobius transformation 
A : Dj nt — > D; nt , where A defined as 

(2) A{z) = Z±±, H 2 -|6| 2 = l. 

bz + a 

This subgroup of Mobius group of selfmaps of the circle is denoted PSL2(M). 

The map / cx t is chosen uniquely via the following normalization: we choose a unique Mobius 
map A, such that / ext o A maps oo to oo, and that its differential carries the real positive axis of 
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the D-plane at infinity to the real positive axis of the T-plane at infinity. Thus the ambiguity in the 
choice of / ex t is eliminated for every F. 

The goal of this construction is to define the map ip which is called the 'fingerprint' (in the 
Teichmuller theory this is known as a 'welding map') of the shape 

(3) V = / ta t /ext e PSL 2 (K)\Diff 

Note, that fc^tiS 1 ) = T, f^ t (T) = S 1 . The fingerprint ip : S 1 — >• S 1 is a real- valued orientation- 
preserving diffeomorphism, and it uniquely identifies the shape T (modulo scaling and rigid trans- 
lations). Due to the Mobius transformation ambiguity in the choice of /i n t, we see by construction 
that ip is a member of the right coset space PSL 2 (IR)\Diff (S 1 ). An example of a shape along with 
four realizations of its fingerprint is given in Figure [2j 

The inverse map from diffeomorphisms to shapes is defined as follows: starting with ip, construct 
an abstract Riemann surface by 'welding' the boundaries of Dj nt and D cxt via ip. The resulting 
Riemann surface must be conformally equivalent to the Riemann sphere. Choose a conformal map 
/ from the welded surface to the sphere taking oo E D ext to itself and having real positive derivative 
there. Let T — /(S 1 ) (for details and the numerical implementation see [30] ) . 

One can equally well define the fingerprint to be 

= /ext' ° /int € Diff(5 1 )/PSL 2 (M), 

which is simply the inverse of our fingerprint. This alternate version is the definition used in [30j . 
However, in this paper we choose right cosets and put the Mobius ambiguity on the left. 

2.2. Weil-Petersson Norm on the Lie Algebra of Diff (S 1 ). The Lie algebra of the group 
Diff(S' 1 ) is given by the vector space Vec(S' 1 ) of smooth periodic vector fields v(9)d/d8 on the 
circle. In [28] it has been shown that the embedding PSL 2 (R)\Diff (S 1 ) <—> T(l) is holomorphic and 
the pullback of the Weil-Petersson metric for v E Vec(S' 1 ) can be expressed as: 

IMIwp = £l» 3 -n|K| 2 

= / Lv(9)v(9)d9. 
Js 1 

Here v(9) = YlnL-oa v ne m6 (where — u_„ for the vector field to be real), and Z = Z\{n — 0, ±1}. 
The Weil-Petersson operator L is an integro-differential operator and it has the form 

(4) L = -H(d 3 e +de). 

Above, H is the periodic Hilbert transform, defined as a convolution with i cot (0/2). 

The null space of the L operator is given by the vector fields whose only Fourier coefficients are 
V-i, vq and v%, i.e. vector fields of the type (a + 6 cos 9 + c sin 9)8/ 89. These vector fields are exactly 
in the Lie algebra sfeO^) of the Lie group PSL20R.). 

2.3. Extending the WP Metric to PSL 2 (R)\Diff (S 1 ). Consider any Lie group G, a subgroup 
H, and let g and f) be their corresponding Lie algebras. 

Any norm || • || on the Lie algebra of G which is zero on the Lie subalgebra of H and which satisfies 
||Ad^(?;)|| = ||u|| for all h e H induces a Riemannian metric on coset space H\G which is invariant 
by all right multiplication maps R g : H\G — > H\G, g E G [2"T] . 

In particular this applies to G = Diff (S 1 ), H — PSL 2 (M) and the above WP norm on vector fields, 
hence it gives the right-invariant WP- Riemannian metric on the coset space PSL 2 (R)\Diff (S 1 ). 
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Consider any two diffeomorphisms i^o^i € PSLa(R)\Diff (S 1 ). The Riemannian distance in- 
duced by the WP norm on vector fields is given by 

(5) L = / \\v(s)\\ WP ds, 

Jo 

where v is a vector field that carries ipo to ip\: 

(6) M0)=M0)+ I v(6,s)ds 

Jo 

This above notation extends for the remainder of this paper: ip\{9) and ipo(@) are the initial and 
target welds (shapes), and the path <fit(0) — 4>{9,t) is a geodesic flow. Vector fields v that minimize 
the distance ^ are geodesies on PSL2(lR)\Diff (S 1 ), and it is a standard fact from variational 
calculus that vector fields v corresponding to geodesies satisfy ||v(s, OIIwp = const. 



3. The geodesic equation 



The Euler-Poincare equation for diffeomorphisms (hereafter 'EPDiff ') is a variant of Euler's equa- 
tions for fluid flow. It describes geodesies on the Lie group of diffeomorphisms of K™ in any right 
invariant metric given on vector fields by \\v\\ 2 = J Rn (Lv,v)dx for some positive definite self-adjoint 
operator L (where (■, ■) is the canonical L 2 pairing). The general EPDiff(]R n ) is derived in [3] and 
has the form 

Q 

—Lv + (v ■ V)(Lv) + divv Lv + Dv l ■ Lv = 0, 

where v is a smooth vector field in R™, V = ^g§^, ■ ■ ■ , gf - ) is the divergence operator, L is a 
self-adjoint differential operator and Dv is a Jacobian matrix. 

The space PSL2(R)\Diff (S" 1 ) that interests us is not a group and is instead a homogeneous space, 
but it has been shown in [BO [5T] that Arnold's formula for geodesies on Lie groups extends to the 
case of a homogeneous spaces H\G. 

Given a path <p t (d) = <p(6,t) in DifT(5 1 ), let v(d,t) = ^(^(fl,*),*) be the scalar vector field 
it defines on a circle and let L be the Weil-Petersson differential operator L — —H(dg + dg). Then 
EPDiff takes the form 

(7) (Lv) t +v.(Lv) e + 2v e .Lv = 0. 

Above, v(8,t) is called the velocity of the path, m(9,t) — Lv(6,t) is the momentum, and this 
equation is the same as introduced in ([!]). We note in particular that the momentum can be a 
distribution. The velocity field v(8,t) lies in the space Vec(S ,1 )/st2(]R), where Vec(S' 1 ) is the space 
of smooth vector fields on the circle and s^M) is the Lie algebra of the group PSL 2 (M). However, the 
momentum m(8) £ Hor, where Hor = {v = ^v^e 1 ^ : iik = 0, k = 0,±1}. Hor is the orthogonal 
complement to s^M) in Vec(5 1 ). We will refer to it as the horizontal space. 

The energy is given by 

E(t) = (m{;t),v(;t)) = \\v(;t)\\wP, 

where (•, •) is the L 2 inner product on [0,27r]. E(t) is constant in time under geodesic flow. The 
v — > m map may be inverted by the relation v(6,t) = G * m(9,t), where G is the Green's function 
G(9) of the WP operator L. The Green's function G(9) is obtained as a solution to LG = Proj(5o), 
where <5o is the Dirac measure centered at 6 = 0, and Proj(^o) is the projection of 5o onto the 
horizontal space Hor. 

The Green's function for the WP operator Q is defined up to addition of any element of s^K), 
a + bsin9 + ccos9 for some constants a,b,c. We normalize G(9) as in [5D] so that it lies in the 
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horizontal space Hor: 

(8) G{9) = (1 - cos 9) log [2(1 - cos 9)] + ? cos 9 



4. Teichons, Singular Solutions of EPDiff 

The EPDiff equation |7| admits momenta solutions that, once initialized as a sum of N Dirac 
measures, remain a sum of N Dirac measures for all time [3[TS]. In reference to this self-similarity 
property, these singular solutions are named teichons (or an iV-teichon). In this paper, the velocity 
field that defines the geodesic between two shapes will be approximated by an iV-teichon. 

For a solution to EPDiff (m), we employ the iV-teichon ansatz 

N 

(9a) m(9,t)=J2Pj(m0-Qj(t)), 

3=1 
N 

(9b) v(9,t) = J2Pj(t)G(e-q j m 



where S = Sq is the origin-centered Dirac mass. Plugging these expressions into EPDiff ([7|), we 
obtain a system of ODEs describing the evolution of the momentum coefficients p k and the teichon 
locations q k : 



(10) 



Pk = ~Pk EjLi P]G'{q k -qj), 
<ik = J2f=iPj G (<lk - qj)- 



As mentioned in Section[3]the momenta m lie in the horizontal space, i.e. m(9, t) must have vanishing 
0th and ±lst Fourier coefficients. Using (9a), we obtain a set of three constraints for (qk,Pk), linear 
in p k : 

N N N 

(11) = 5>e^ = 5>e-^ = 0. 

If they are satisfied at time t — they will be satisfied for all t. The teichons never collide: i.e. 
the teichon locations q k retain their initial ordering on S 1 for all time. However, it is known that 
most initial configurations lead to an exponential decay of teichon separation, and an exponential 



increase in the momentum coefficients |20j . This introduces numerical difficulties in solving (10) 



5. Numerical method 

Our numerical method has many components, and each of them requires some discussion. We 
proceed through these components as follows: we discuss construction of welding maps in Section [5.1| 
The overarching shooting method is presented in Section [5T2| and the cross-ratio matching term is 



introduced in Section 5.3 The construction of this matching term is nontrivial, involving a Delaunay 



triangulation of points in the complex plane, and Section 5.4 highlights these considerations. The 
(standard Euclidean) gradient of the matching term is necessary in order to update our initial 
teichon guess, and computation of the gradient is discussed in Section |5.5| It is well-known that 
the gradient is not the optimal search direction for optimization on Riemannian manifolds, and in 
addition the computed gradient does not satisfy the momentum admissibility conditions ( flT] ); in 



Section 5.6 we transform and project the gradient to address these concerns. Finally, Section 5.7 
discusses refinement strategies for obtaining good initial guesses, and Algorithm Listing [l] presents 
the full algorithm. 
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5.1. Generating fingerprints. We first consider the task of constructing fingerprints. In practice, 
motivated especially by computer vision, we will be given an ordered collection of points {z m }^ =1 C 
C lying on some closed curve in the complex plane. (This is our shape.) Naturally this does not 
define a continuous curve in the plane, but this situation is realistic in an application setting. 

From this discrete data, we use conformal welding to construct a continuous weld ifi € PSL2 (R)\Diff (S 1 ) 
This task requires only the ability to construct conformal maps between the unit disc and a region in 
the complex plane defined by the z m . While many methods are suitable (in particular the methods 
described in [30]) we choose the Zipper algorithm [24] , which computes the map via a discretization 
of the Loewner differential equation. One main strength of the algorithm is its sequential nature 
- computation of the entire map is a method that simply iterates over the index m on the input 
points z m in an explicit way. In particular, the total work required is 0(M 2 ). 

This should be contrasted with a popular competitor, numerical Schwarz-Christoffel mapping [9, 8 
requiring the solution of a size-M optimization problem, which can exhibit convergence problems, 
which are exacerbated when the fingerprint has very large, or very small derivatives. The sequential 
nature of Zipper means that one computes the full conformal map as a composition of intermediate 
maps; when the fingerprint has exponentially large derivatives, such a compositional strategy is more 
computationally robust. Our experience indicates that Zipper algorithm is more efficient, and more 
resilient, for generating welding maps. We refer the reader to [24j for details on implementation of 
the Zipper algorithm, where convergence in the Hausdorff metric is proven. 

The situation of very large, or very small derivative values for a fingerprint is called "crowding" 
in the literature: Recall that conformal welds given by ^ are a composition of conformal maps /i nt 
and / ox t . Given a conformal map / whose domain is the unit disk Dj nt , it is crowded if 

= max xeS i \f(x)\ 
minxes 1 I/' 0*01 

is large, where the definition of 'large' depends on the finite-precision arithmetic being performed. A 
good rule-of-thumb for 'large' is when R is inversely proportional to machine precision, R ~ e^Ich- 
A reasonable characterization of a crowded weld then is when the product of the R values for the 
interior and exterior maps is on the order e~^ ch . Conformal welds ijj are often stored as point- 
evaluations (O cx t,mi $int,m) that are computed from boundary values of the conformal maps /j nt and 
/ oxt . Any point-evaluations that are sampled in regions of S 1 where either map is crowded will 
coalesce to machine precision; that is, 

^int,m+l ^int,m ^cxt,m+l ^ext,m 

j ^raachi Or - j ~ ^mach- 



When this happens, the weld ip is effectively not a diffeomorphism to machine precision, and it is 
difficult to accurately compute particle locations in crowded regimes. In the context of flow under 
EPDiff, particles 9 m may, at t = 0, start in an uncrowded configuration (that is, the t = weld is 
not crowded), and then flow to a crowded configuration (that is, the t = 1 weld is crowded). 

Our algorithm does not ameliorate the underlying problem with crowding; indeed both Schwartz- 
Christoffel mapping and the Zipper algorithm do not produce accurate results for crowded shapes^ 
However, we mention again that our experience is that Zipper is more robust when the shape is 
crowded. 

Finally, we remark that although the input to the Zipper algorithm is a discrete set of points, 
the output is a continuous welding map %p. Therefore in the sequel we continue to speak about 
continuous welding maps. 



^We acknowledge the possibility that a solution is given in 1101 . but an application of this method to conformal welding 
is a separate, self-contained project in itself. 
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5.2. Shooting method. We want to compute the geodesic between the two shapes, given by finger- 
prints ipQ and ipi . In other words we seek to find the velocity field v corresponding to a geodesic for 
([7]) such that the diffeomorphic evolution defined by Q satisfies the prescribed boundary conditions 
4>(t = 0) = ipQ and <f>{t = 1) = ipi. We will solve this problem by approximating the velocity field by 
an evolving iV-teichon solution to EPDiff. 

Since |7]) specifies the evolution of an initial velocity field, the goal then is to find the initial 
positions of teichons, <7fc(0), and initial teichon strengths, Pfc(0), such that the resulting velocity field 
will carry a template fingerprint %pQ as near as possible to the target ipx. The task of determining 
initial data to satisfy a two-point boundary value problem is well-studied and one of the more popular 
numerical methods to compute a solution is the shooting method [29 . Diffeomorphic matching in 
the context of shapes has also seen the recent application of shooting methods [55] . 

The idea of the shooting method is the following: start with an initial guess for the 5fc(0),pfc(0), 
construct the initial momentum m((9, 0) = ^2pk(0)G(9 — g^(0)), and then solve forward the equation 



(10 1 to obtain time evolution of qk(t)iPk{t)- This in turn will produce a time- varying velocity field 



v{9 1 t), which we integrate via the equation 

(12a) v(9,t) = ^(^ 1 (9,t),t), 

(12b) ^,t = 0)=Vo- 

to obtain (f>(9, t = 1). The final computed fingerprint <p(9, t = 1), is compared with the target finger- 
print, -01. Based on this comparison, we modify the initial shot configuration {pfc(0), 9fc(0)}^ =1 and 
repeat the process. Because the pk must satisfy constraints that depend on q^, our shooting method 
only changes the teichon momenta pk] changing is certainly possible but requires admissibility 
constraints whose application is more involved. We have found that varying only the momenta 
allows us to represent a large variety of shapes. 

In practice, we cannot match ipi up to infinite precision, so we resort to inexact matching via 
some discrete set of landmark points. Given the discussion from Section pTT| it is sensible to choose 
M landmarks corresponding to the images of the terminal shape samples z m . (Here, terminal means 
the target shape at t = 1.) We recall that fingerprints ip are defined through conformal maps /: 
"00 = (/o~int /o.exOls 1 ! and similarly for ipx. We consider landmark points 9 m g [0, 2tt) on the 
exterior defined by 

(13) exp(i6» m ) = /^ st (z m ). 

We track these landmarks as they flow from the initial shape: let a m = ipo(9 m ) — argo/~ ; ^ t o 
fo,cxt{z m )- We flow these landmarks using ( 12a[ ) to t = 1 and in principle we wish to compare their 



locations with the exact terminal locations ipi(0 m ) = j x l n ti z m)- The t — 1 images of a m under an 
A^-teichon evolution are determined by 



N 

p n G(a m -q n ), 



n=l 



(14) 

The particular choice that a m (0) = ipo{9 m ) is not the only choice one could make, and we do not 
claim it is optimal; however, our results indicate that such a choice performs quite well in many 
situations. The choice of landmark locations a m and the t = teichon configuration q^ need not 
be related. With these M landmark locations a m (l), we must compare fingerprints at t = 1. We 
thus need a matching functional E2{ipi 1 ip{0 1 1)) that reflects the closeness of the fingerprints. The 
relation between the terminal shape samples z m and the landmarks 9 m is shown in Figure [l] 



5.3. Matching fingerprints. The shooting method relies on computation of a matching term, 
which we call E. This matching function compares the fingerprint computed with an iV-tcichon 
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Samples on S 1 (exterior) 



arg(-) 



2- 




Samples on S 1 (interior) 



arg(-) 



Figure 1. For a terminal shape, relationship between the desired landmark positions ipi(8 m ), 
and the welding map generated from the shape samples z m . The initial landmark locations a m = 
i>o(fim) are not shown. 



evolution with a target fingerprint. Note that the standard type of matching functional for geodesic 
shooting in this context would take the form 



(15) 



E= / \\v(;t)\\^ p dt + Xj2 WlW-OmW) 
J m= 1 



= Ei + i?2- 



for the landmark choices {O m }m=i C S , and a scalar A > 0. The weight A defines the relative 
importance between the first energy term, and the second matching term. A shooting procedure 
would iterate on the initial data in an attempt to minimize this functional. Our version employs two 
variations. First, we are not minimizing energy: geodesies on PSL2(R)\Diff (S 1 ) are unique [13] so 
that the value of the energy is irrelevant since only one path exists between ip and V'l- Therefore 



we entirely omit the first energy term E\ in (151 



Our method also employs a different landmark matching term. We cannot directly employ an 
^ 2 -type distance as given in ( 15 ) because of the Mobius invariance of welding maps. For comparison, 



we show four different welding maps associated with the same shape in Figure [2j We want any 
matching term that we devise to assign zero distance between any pair of welding maps in the 
figure. However, using an £ 2 type distance to compare them is clearly misleading. In particular, 
given any fingerprint ipi and e > 0, there is a Mobius map A m e PSX2OR) such that the pointwise 
distance between and A m o ipi at 6 m is within e of the maximum matching distance: 

(M0 m ) - {A m o Vi) {e m )f > 4tt 2 - e. 

Given ip and ipi, we cannot determine the appropriate Mobius maps so that ipo an d ipi have 
the "same" normalization until we have already computed the geodesic connecting them. Since the 
matching term should be insensitive to self-maps of the disk, we amend the landmark matching term 
to be robust with respect to projective transformations on C. The most general projective invariant 
quantity of a 4-tuple of points in the complex plane is the cross-ratio and forms the basis for our 
matching term. 
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Figure 2. Left: an ellipse. Right: various welding map representatives of the ellipse shape 
corresponding to the same equivalence class in PSL2 (M)\Diff (S 1 ). 



Let (zi, z 2l Z3, Z4) C C be a 4-tuple of points. The cross-ratio of these four points is defined by 

(16) C(zi,z 2 , z 3 ,z 4 ) = r. 

(z 2 - z 3 )(zi - z 4 ) 

C is invariant under the Mobius transformations: C(z\, z 2 , 23, z 4 ) — C(A(zi), A{z 2 ) 1 A(z^), A(z 4 )) 
for any Mobius map A. One can also show that if Zj € S 1 for all j, then Cel. 

The invariance of the cross-ratio under Mobius transformations will allow us to compare finger- 
prints that have different Mobius normalizations. After evolving according to EPDiff, we have M 
landmark locations on S 1 that specify pre-images of shape vertices under tpi nt . This suggests that we 
can only resolve the shape up to these M vertices, and that furthermore we can only use cross-ratios 
of these pre-images in order to have the Mobius invariance. 

A method to uniquely encode information about a polygon with M vertices has been proposed 
in [10]. The basic idea is that a polygon is uniquely identifiable if the vertex angles of a poly- 
gon are specified in conjunction with M — 3 carefully chosen cross-ratios of quadrilaterals. These 
quadrilaterals are constructed from the Delaunay triangulation of the polygon. 



5.4. Delaunay triangulation. Let P be a simple polygon. A triangulation of P is a division of 
P into non-degenerate triangles whose vertices are vertices of P. The triangles intersect only at a 
vertex or at an entire edge. A Delaunay triangulation of polygon P is a triangulation such that no 
point in P is inside the circumcircle of any triangle in the Delaunay triangulation. Triangle edges 
of the triangulation that are not polygon edges are called diagonals. 

It is known [5] that every P has at least one Delaunay triangulation with the following property. 
If d is the diagonal, let Q(d) be the quadrilateral, composed of the union of two triangles on either 
side of d. Then the sum of two opposite interior angles of Q(d) that are split by d is at least n. A 
Delaunay triangulation of an n-polygon P can be computed in 0(n 2 ) steps. In our implementation 
we have used the MATLAB function delaunay. 

It is well known that any n-vertex simple polygon P has a triangulation consisting of n — 2 
triangles. In addition it has exactly n — 3 distinct diagonals. The vertices of a quadrilateral Q(d) 
associated with each diagonal are used in the computation of the cross-ratios, see Figure |4j More 
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specifically, it is shown in [TDJ that using cross-ratios computed using this choice of quadrilaterals 
uniquely characterizes the original polygon P. 

With this in mind, we let {ik,i, ik,3, ik,4}k=i denote the K 4-tuples of cross-ratios from [ID] . 
If we have M points for our shape, i k j is an index taking a value in 1, . . . , M indicating which point 
to use in the fcth cross-ratio. We take the matching functional to be the relative error of the discrete 
£ 2 cross-ratio difference: 



(17) [Coeo«,l)] W „,, 



Ik, 2 1 "»fc,3 ' "lk,i i 



K h^\ ^ ° e0 ^ 6 ^,n e ik,2^lk,n e ikA) 



where e(x) = exp(ia;) is the complex exponential, and the 8 m are the points defined by (13). 



We emphasize that the construction (17) for the matching term is automated: the identification 
of the required quadrilaterals making up the cross-ratio term is automatically computed using the 
Delaunay triangulation of the shape. After a one-time run of the Delaunay triangulation, the formula 
(17) is also explicit: the indices ik j are known and stored. 



5.5. Gradient with respect to p(0). In order to adjust the initial momenta to reach the target, 
we need to compute the gradient of the matching functional, E, with respect to p{0). Let po denote 
p(0), then the gradient of the energy is given by 



(18) 



BE _ dE da 
dpo da dpo ' 



where a is the length- M vector of landmarks solving (14) at time t — 1. Define (3 m = 7r n — 
and Xn — i each of which is a 1 x N vector for m 



1. 



, M and n = 1 , 



A dp™ 
~ dpo 

N. These parameters 



may be computed by a system derived from ( 10 ) and ( 14 ) 

TV 

/?m [ir k G(a m - q k 



(19a) 
(19b) 
(19c) 



d/3„ 



dt 

dir n 
dt 

dXn 
dt 



PkG'(a m - q k )(/3 m - Xk)} 

N 

^n^PkG'(q n - q k ) - p n ^ [KkG'(q n - qk) +PkG"(q n - qk)(Xr. 

k=l k=l 

^2 \' K kG{q n - q k ) +p n G'{q n - q k )(Xn - Xk)} , 



k=l 
N 



Xk) 



N 



k=l 



where we assign G"(0) = 0. The full system (19), (14), and (10) can be solved in parallel to determine 
/3 m (l) = to be used in (18). One can explicitly computel^ from (16), and the definition of the 
complex exponential e(-). Thus the energy gradient (18) is computable. 



5.6. Optimization with the gradient. There are three tasks yet to be accomplished before we 
can update the initial momentum distribution: 

• dE/dpo is not an admissible momentum distribution, so we must project it into the appro- 
priate space 

• on non-Euclidean Riemannian manifolds the gradient does not point in the direction of 
steepest ascent; we require the natural gradient 

• gradient descent is the most basic of optimization methods; we employ a nonlinear conjugate 
gradient update to accelerate convergence 

This subsection discusses these considerations. 
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5.6.1. Projecting the gradient. For the Teichon evolution system ( 10 1 to be valid one needs to have a 



bijection between the Lie algebra g and its dual Q* . This bijection is provided by the Weil-Petersson 
operator L and it's inverse, convolution with the Green's function. These operators are a bijection 
only on the horizontal space Hor = {w € Vec(S' 1 ) : Wk = 0, k = 0,±1 and ||u>||wp < oo}, where 
the Wk are Fourier coefficients of the periodic function w{9) defined on the circle. In other words: 

L : Hor -> Hor* , 

G : Hor* Hor 

LG = G * L8 = Proj(S). 

Here, Proj((5) is the projection of the delta function onto the horizontal subspace Hor. Therefore 
any updates we perform to the initial momentum distribution must happen on this space. In order 
for a momentum field m £ Vec(S' 1 ) to lie in Hor*, it must likewise have vanishing 0,±1 Fourier 
coefficients. This defines the three constraints given by (11) as discussed in Section [lj 



While our starting guess for p(0) will satisfy the conditions (111, there is no guarantee that the 



update vector from (|24|) will satisfy those constraints. We must therefore obtain an clement 



from Hor' 



Let Ap" e 



BE 

S lven 9^' 

8E 
dpo ' 



which does not represent a member of Hor*. We project this update vector 



into the space admissible updates: those that satisfy ( 11 ). We proceed by computing the VFP-closest 
member of Hor* to Ap 71 



Let the 3 x N matrix F define the admissibility constraints: 

F{q(0)) = ( cosgi(O) ... cosgAr(O) 
\singi(0) ... singAr(O) 



We wish to find an update vector Ap satisfying 

Ap 



(20) 
where 



argmin 1 1 Ap 

F(q(0))Ap=0 



Ap n 



\WP* 



\wp* is the norm on Hor* induced by the norm || • || 

\\Ap\\ WP * = \\G * Ap\\ WP . 



rp on Hor: 



Let G(g(0)) be the NxN Gram matrix for {Proj((5 9ji ( ))}„ =1 , whose entries are given by Gi t j(q(0)) = 
G(qi(0) — f? j ( ) ) , where G(-) is the Green's function wh. Then (20) can be written as minimization 



of a quadratic objective subject to a linear constraint. The solution is a vector Ap given by 



(21) 

where above G :- 



Ap ■ 

G(q(0)) and F : 



I -G- 1 F T (FG- 1 F T )- 1 F 
F(q(0)). 



Ap, 



5.6.2. The natural gradient. It is well-known that performing gradient descent on Riemannian man- 
ifolds with the standard gradient is not the optimal gradient update. It is much more effective to 
use the natural gradient (see, e.g., [2]). The idea behind the natural gradient is the following: given 
an A^-teichon p on the space Hor*, suppose Ap is the standard gradient direction (satisfying the 
admissibility constraints from the previous section) . We wish to move in the direction that decreases 
the objective most, and so we must solve the problem 

Ap = argmax (p, Ap)wp* = argmax p T GAp 

IIpIIwp*= 1 IIpIIwp*=1 

where the constraint ||p||wp* = 1 may be written as p T Gp — 1 where G is the Gram matrix for the 
iV-teichon configuration of p. The solution to this problem is given by 



(22) 



Ap = XGAp, 
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for some normalizing constant A, which we hereafter set to unity. The new gradient direction 



Ap is called the natural gradient. Combining this with (21 1, the full update vector Ap given the 
unconstrained gradient is 

dE 



(23) 



Ap 



G - F T {FG- 1 F T y 1 F 



dp 



In the sequel we will refer to the momentum update (23) that is both admissible (satisfying (111) 



and natural (given by (22)) as the proper gradient 



5.6.3. Updating the shooting direction. With the proper gradient vector Ap given by (23), we can 



proceed with standard optimization methods. The most straightforward is steepest descent: The 
update for the vector p(0) at each iteration is given by 

(24) p(0) <- p(0) - eAp 

where the choice e determines how far along the gradient direction Ap we update. 

Convergence with gradient descent often stagnates when the iterative path taken by solutions 
follows a narrow valley; in such cases more sophisticated methods are required to render the iteration 
computationally efficient. A nonlinear conjugate gradient method is one such alternative. We employ 
the Polak-Ribiere method, which makes use of the update vectors from the previous iteration. Let 



p n denote the proper gradient (23) at iteration n. Compute 

(Pn, Pn — Pn-l}wP* 



f3 = max < 



(25) 



Set the update direction at iteration n to u n 
where e* satisfies 



(Pn-li Pn-l)WP* 

p„ + /3oj„_i. Perform the update p(0) 



P (0) 



■e u* 



(26) 



argmini?(p(0) — ecu n 

e>0 



To begin the iteration process, the first iteration is performed as a standard gradient descent update. 

We utilize the standard nonlinear optimization tricks: standard gradient descent is performed 
at the initial stages to iterate close to a basin of attraction; then nonlinear conjugate gradient is 
employed to quickly converge. 



5.7. Initial configuration guesses. In all of our experiments, we specify the number of teichons 
N, which is fixed throughout iteration. We also determine the initial teichon configuration q(0) by 
spacing the N teichons equidistantly on S 1 . This is not the optimal choice, especially for shapes 
with important features concentrated in a particular location, but our tests indicate that this is not 
a bad choice for many non-crowded shapes. 

In taking the initial choice for p(0) to be the zero vector, we have found that convergence takes 
an inordinate amount of time, or the iteration stagnates and convergence is not observed at all. One 
well-known way to combat this is to choose a better initial guess. We do this by first solving the 
minimization problem on a subset of the full collection of cross-ratios. We identify this subset by 
choosing cross-ratios that resolve coarse features of the shape. 



We therefore implement a coarse-to-fine approach in the matching functional (17). We compute 
the Delaunay triangulation on a coarse subset of the vertices of a given polygon P. For example, if a 
polygon P has M — 128 vertices, we compute a triangulation on dyadic subsets consisting of 8, 16, 
32 and 64 points, and finally the full set of 128 points. We use a zero initial guess for p(0) and run 
the minimization algorithm for the cross-ratios identified by the 8-point Delaunay Triangulation. 
The solution of the 8-point problem is used as the initial guess for the 16-point problem, and so 
on. Assuming a nested choice of points (i.e. the 8 points are a proper subset of the 16 points, etc.) 
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then we progressively build up the choice of cross-ratios until we utilize those for the full set of 128 
points. 





— ■'. ... 



.A 



Figure 3. Coarse to fine Delaunay triangulation used in the matching functional \17\ for the 
shooting method. Triangulation is depicted on the subsets of 8, 16, 32, 64 points and the full set of 
128 points of the fish. 




Figure 4. The 4-tuples used in computing of cross ratios in the matching functional |[17l arc 
connected by solid lines. Dotted lines depict Delaunay triangulation on the coarsest subset of N 
points representing the fish (N = 8). For N points used in computing the cross ratios we get N — 3 
four-tuples. 



We show the Delaunay triangulation for a fish shape at various refinement stages in Figure [3] 
Our choice of using a dyadic refinement strategy is not the only possibility, and we have chosen it 
mainly for convenience. In our results we employ between 4 and 5 refinement stages. 
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We summarize the full shooting procedure outlined in this section in Algorithm [I] 

Algorithm 1 Shooting algorithm for minimization of the cross-ratio objective ( p~Tj ) . 

Input: number of teichons N and M ordered shape samples z m € C 
Determine number of refinement stages S: Section [5. 7| 
for Stage s — 1, ... S do 



Compute Delaunay triangulation and cross-ratios for 2^ s S ^M shape samples: Section 5.4 



Set initial guess for p(0) as solution to previous stage: Section 5.7 
while not converged do 



Compute objective E% (17) and the standard gradient (18) using (10), (14), (16), and (19): 
Sections [531 [531 



Compute proper gradient using (23): Sections 5.6.1 5.6.2 



Perform either steepest descent (24), or conjugate gradient descent (|25|), (|26|: Section 5.6.3 
end while 
end for 

Output: A-teichon configuration q(0), p(0) 



6. Numerical Results 



In this section we present various numerical results that demonstrate the efficacy of our method. 
For all evolutions we employ A = 100 teichons equispaced at t = 0. Unless noted otherwise we 
iterate until the value of the objective (17) is no greater than 10~ 4 , and frequently is 0(1O -6 ) at 
convergence. The number M of landmarks a m varies with the shape data, but usually takes the 
value M = 128. 



6.1. Aspect ratio for an ellipse. We first investigate the distance on PSLi2(R)\Diff (S 1 ) between 
a circle and an ellipse of certain aspect ratio. With M = 100 matching points, Figure [5] graphs 
the distance on T(l) between a circle and ellipses of aspect ratios from 1 to 6; these results are 
visually indistinguishable from those found with an energy minimization algorithm in |11| , providing 
supporting evidence for the accuracy of the algorithm. The figure also suggests that the asymptotic 
ratio between the geodesic length and the aspect ratio of the ellipse is linear and the slope is 
approximately 0.69. We note however that data for larger aspect ratios is necessary in order to verify 
this result. Unfortunately, numerical crowding prevents us from accurately computing geodesies for 
higher aspect ratios. 

6.2. Hyperbolicity test. As has been mentioned in Section [T] all sectional curvatures of the Weil- 
Petersson metric are negative. In order to verify this numerically, we verify that the angle sum of a 
triangle on T(l) with the WP metric is less than ir. We choose the three vertices of each triangle to 
be rotated ellipses of fixed aspect ratio. For a template ellipse with semimajor axis aligned with the 
horizontal axis, we choose the three vertices to be ellipses of rotations 0, 27r/3, and 47r/3; we label 
these vertices s%, S2, and S3, respectively. We compute geodesies with M = 128 matching points 
between two vertices of the triangle, and compute angles between these geodesies. 

Let Vij £ Hor denote the velocity field that pushes vertex Sj to vertex Sj. The angles at the 
three vertices of this triangle are computed according to the formula: 

\Vi,i@l, ViAQl/WP 

Ql^ — ' ^ 

H w i,i©l||wp|| w i,iel||wP 

where © and G denotes modular addition/subtraction on the set {1, 2, 3}. We then compute X^=i a i> 
which is the angle sum of the triangle. The graph of this sum versus the aspect ratio of template 
ellipse is depicted in Figure [6] The aspect ratio of ellipses varied from 1 to 2.2. 
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1 2 3 4 5 6 



Aspect ratio 

Figure 5. Left: WP distance from circle to an ellipse vs. aspect ratio of the ellipse computed 
using the shooting algorithm presented in this paper (solid line) and using an energy minimization 
algorithm from (dashed line). (The lines overlap.) Right: Snapshots of the computed geodesic 
for aspect ratio 3 at equidistant points along the path. 

Once again we compare our results against those computed from an entirely different algorithm 
in This comparison is shown in Figure |6j and in contrast to the previous test, there is now a 
noticeable difference in the results: The values differ by about 2%. We can attribute this difference 
to many factors. First, the sample points on the shape that we use in this algorithm are not the 
same points that are used in [TT]; this results in small differences in the conformal welds used for 
matching. Second, our algorithm uses M = 128 matching points, whereas the energy minimization 
algorithm from |llj used 150 points. Finally, our shooting method produces a numerically exact 
geodesic, but does not match the endpoint condition exactly; the energy minimization method from 
[llj produces an approximate geodesic that matches the endpoint condition exactly to numerical 
precision. Thus, the small differences shown in the results are not altogether surprising. 

6.3. A hippocampus slice. We consider flowing from a circle to a planar slice of the human 
hippocampus in Figure [7] given by M = 128 sample points. (For details, refer to [22].) We shoot 
with 100 equidistant teichons and at algorithm termination the objective function is less than 10~ 4 . 
Figure [7] shows that the shooting matches the target very well. Figure [8] displays the shape evolution 
along with a contour plot of the velocity field that pushes the circle to the hippocampus slice. We 
also track the evolution of 5 landmarks on the shape. The hippocampus slice is a relatively easy 
shape: the fingerprint does not exhibit crowding and so our computation of the geodesic flow (and 
the gradient) is accurate. 

6.4. Flow between shapes. Our method does not rely on the initial shape being circular - it 
is likewise possible to flow between non-circular shapes with no change to the algorithm. The 
MPEG- 7 CE-shape-1 collection of planar shapes [I] is a database of shapes commonly used in 
classification routines. Our immediate goal here is not classification, but to illustrate the applicability 
of our algorithm to realistic shapes, we choose two shapes with non-crowded welding maps from this 
database and show the Teichmuller evolution between them obtained from the shooting algorithm 
with M = 128 matching points, see Figure [9] 

6.5. A fish contour. We finally consider a more complex shape: an outline of the fish given in 
Figure [10] from M = 128 sample points. The initial shape is a circle, and the initial teichon 
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Aspect ratio 



Figure 6. Left: Angle sums for triangles on T(l) computed using the shooting algorithm of 
this paper (solid line) and the minimization algorithm from 1111 (dotted line). The three vertices 
for each triangle are formed by rotating a template ellipse of given aspect ratio. Right: sample 
evolution for ellipse of aspect ratio 1.6. An artificial shift is employed to make the evolution clearer. 



Figure 7. Terminal shape with shooting (solid line) versus landmarks used for the cross-ratio 
objective (solid dots). 



configuration is given by 100 equidistant teichons. The result of the matching is given in Figure [TO] 
As with the hippocampus slices, we show the resulting WP geodesic from a circle to the fish on the 



right panel of Figure 11 and a contour plot of the velocity is shown on the left panel. The objective 
value at termination of the algorithm is 5 x 10~ 4 . We see immediately that there are limitations to 
this algorithm: the welding map for the fish outline suffers from severe crowding. 

For the particular chart we have chosen for the fingerprint ip\, the landmarks on the tail of the 
fish are separated by a distance of 0(1O~ 8 ). This crowding of the shape landmarks a m implies a 
similar crowding of the teichon positions qk that push the landmarks. When teichon positions qk 
and qk+i are 0(1O~ 8 ), we can no longer accurately integrate EPDiff. To understand why, consider 
first the Gram matrix G for teichon positions qk (see (21)), which has entries Gk,i = G(qk — qi) 1 
with G(-) being the WP Green's function (J8|. From the explicit form of the Green's function, one 
can show that for small arguments A9, 



G(Ad) 



A9 2 



log(A<9 2 
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Time t 

Figure 8. Left: Contour plot for velocity field v(9, t) evolving a circle (t = 0) to a hippocampus 
slice (t = 1). Right: The resulting shape evolution. Shape snapshots are shown for t = n/20 for 
n = 0, . . . , 20. An artificial scaling is employed to make the evolution clearer. The path length on 
Tciclimullcr space is 1.968. 




Figure 9. Shooting from a bell to a blob, two shapes in the MPEG-7 CE-Shape-1 database. 
Left: Terminal shape (blob) with shooting (solid line) versus landmarks used for the cross-ratio 
objective (solid dots). Right: Shape evolution with landmark flow (solid dots). An artificial scaling 
is employed to make the evolution clearer. The path length on Teichmuller space is 6.362. 

Then when qk — qh+i is 0(1O -8 ), we have Gk,k+i = G{qu — Qk+i) ~ \ + 0(1CT 16 ), and in double- 
precision arithmetic where we have implemented this code, floating-point truncation error causes 
this matrix entry to coincide with Gk.k = G(0) = \. This means that the Gram matrix is singular 
to numerical precision. Thus we cannot accurately evaluate the right-hand side of EPDiff given by 
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Figure 10. Terminal shape with shooting (solid line) versus landmarks used for the cross-ratio 
objective (solid dots). 




Geodesic flow to a fish 



^ C 3 ° o a a a - 



0.4 0.6 
Time t 




Figure 11. Left: Contour plot for velocity field v(9, t) evolving a circle (t = 0) to a fish (t = 1). 
Right: The resulting shape evolution. Shape snapshots are shown for t = n/20 for n = 0, . . . , 20. 
An artificial scaling is employed to make the evolution clearer. The path length on Teichmiiller 
space is 7.152. 



(10 1 and (14 1, and also cannot accurately compute the gradient using the EPDiff-derived system 



( 19 ). Therefore the algorithm begins to break down at this point in the sense that we cannot resolve 



features that require teichons to flow so close to one another on S 1 . 

In general, when teichons flow very close to one another (relative to machine precision), we observe 
signature failures of the algorithm due to finite precision through various diagnostics: 



• the computed WP norm of the iV-teichon is not constant in time t, or becomes negative, 

• teichon locations qk (or landmarks a m ) cross each other, 

• stepping in the direction of the proper gradient does not decrease the objective. 



The first two issues can normally be ameliorated by decreasing the At time-stepping parameter used 
to integrate (10), but the third issue is usually difficult to resolve in an automated fashion. 
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7. Conclusions 

In this paper we have demonstrated an efficient method for computing geodesies on the coset 
space PSL2(IR)\Dirr(5 1 ), a dense subset of the universal Teichmuller space T(l), with the Weil- 
Petersson metric via shooting. The geodesies are found by approximating the velocity field with 
an ansatz solution, an A-Teichon. The fact that an A^-Teichon solution remains an 7V-Teichon 
under geodesic flow allows us to accurately compute these geodesies. A matching term is employed 
to guide the initial guess; cross-ratios make up the matching term to correctly identify disparities 
between equivalence classes on the coset space PSLi2(R)\Diff (S 1 ). We are able to use a nonlinear 
optimization algorithm to converge to geodesies on this space. However, our method still suffers 
from the well-known crowding phenomenon, which prevents us from computing geodesies between 
shapes with elongated features. Even if a crowded welding map can be accurately computed, the 
geodesic evolution becomes inaccurate when particles flow to within ^/e mac h , where e mac h is machine 
precision for floating-point computations. Nevertheless, there is a wide range of non-crowded shapes 
for which our algorithm is effective. 

To our knowledge this is one of only a few numerical algorithms that can reliably compute the 
Weil-Petersson geodesies on the T(l). It performs much better than the energy minimization method 
proposed in |30) , and is competitive with the recent approach [11] . An additional advantage to the 
proposed method is that it is a shooting method: a proper geodesic is always produced (up to the 
precision of the forward integration scheme). Future work will employ this method for consistent 
comparison of shapes in a database. The uniqueness of geodesies implies that consistency is ensured 
by the unique initial momentum that is assigned to each shape via the tangent space linearization. 
Moreover, this allows us to find unique Karcher mean and perform well-posed statistics on the shape 
space 18J . In [22] we have employed the described method to study the database of hippocampus 
of patients with dementia and healthy controls. 
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